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Abstract: Diffuse reflectance spectroscopy is one of the simplest and 
widely used techniques for the non-invasive study of biological tissues but 
no exact analytical solution exists for the problem of diffuse reflectance 
from turbid media such as biological tissues. In this work, a general 
treatment of the problem of diffuse reflectance from a homogeneous semi- 
infinite turbid medium is presented using Monte Carlo simulations. Based 
on the results of the Monte Carlo method, simple semi-empirical analytical 
solutions are developed valid for a wide range of collection geometries 
corresponding to various optical detector diameters. This approach may be 
useful for the quick and accurate modeling of diffuse reflectance from 
tissues. 
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1. Introduction 

Many studies in biomedical optics have employed diffuse reflectance spectroscopy for the 
non-invasive characterization and analysis of biological tissues. The technique is simple and 
straightforward in its application and can be implemented with relatively inexpensive 
instrumentation. The measured spectra contain information about the absorption and 
scattering properties of tissues under study, which can potentially be used to obtain 
biochemical, morphological, and histochemical information. However, no exact analytical 
solution exists for the problem of light propagation in turbid media such as biological tissues. 
Early attempts to model diffuse reflectance from biological tissues were characterized by 
cumbersome analytical expressions based on the diffusion approximation [1,2]. Recently, 
several studies have introduced relatively simple models, which are designed for specific 
probe geometries [3,4], with some of them being empirical or semi -empirical in nature while 
others are based on a relatively simple solution of the diffusion approximation utilizing the 
method of images [5]. On the other hand, numerical methods such as Monte Carlo (MC) 
simulations [6] are accurate and they are not restricted to any particular light delivery and 
collection geometry but require long computation times and they are difficult to apply in the 
inverse mode since they provide no analytical solution. Many researchers have tried to 
address these shortcomings by coming up with various ways to accelerate the execution time 
of MC simulations [7-9]. More recently, significant acceleration using graphics processing 
units (GPU) has also been possible [10-12]. Overall, these improvements yield an advantage 
of several orders of magnitude faster execution, making MC simulations a very appealing 
technique. At least one study has recently been published utilizing a MC method in the 
forward mode in lieu of an analytical solution [13]. 

In this work, MC simulations are used to solve the problem of diffuse reflectance from a 
homogeneous semi-infinite turbid medium for arbitrary detector diameter geometry. To 
overcome the lack of an analytical solution using the MC method, a simple semi-empirical 
solution is developed based on the results of the MC simulations, which can be used to 
describe the MC results quite well. In addition, several interesting properties of diffuse 
reflectance are discussed and pointed out in terms of the semi-empirical modeling approach 
and the particular detector diameter. The overall modeling approach and results may be useful 
for the further modeling of reflectance from more complicated tissue geometries, including 
multi-layer morphological structures. 
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2. Methods 



2.1. Monte Carlo simulations 

Simulations of light propagation in turbid media such as biological tissues were performed 
using the publicly available MCML code [6]. A one-layer semi -infinite tissue model was 
assumed (implemented as a 1000 mm thick tissue layer) and incident light was modeled as an 
infinitely thin beam perpendicular to tissue. Reflectance was scored up to a distance of 10 mm 
away from the incident beam with a step size of 0.5 mm so as to simulate a wide range of 
delivery/collection separations corresponding to detectors with various diameters. Reflectance 
at various detector radii was calculated by integrating the total reflectance from r = 0 up to the 
actual radius r of the detector. Total reflectance (corresponding to r = oo) was also scored. 
Reflectance was scored at various photon exit angles (with respect to the perpendicular on the 
tissue surface) corresponding to detectors with different numerical apertures. All possible exit 
angles from 0 to 90 degrees were covered in 6 different bins, each one 15 degrees wide. The 
refractive index of tissue was assumed equal to 1.4 and the refractive index of the medium 
above was assumed equal to 1.0. For simplicity, the refractive index of the detector was also 
assumed equal to 1.0 and was otherwise assumed to be an ideal detector (i.e. with no 
absorption, reflection, etc.). The scattering anisotropy coefficient of tissue was assume to be g 
= 0.9 (Heaney-Greenstein scattering phase function [6]). Table 1 summarizes the distinct 
values of the absorption and reduced scattering coefficient used in the simulations. These 
cover the entire typical relevant range for soft biological tissues. There were 5 x 13 = 65 runs 
of the MCML code in total, each one for each pair of absorption and reduced scattering 
coefficient values. Every MCML run was performed with 10 million incident photons. All 
reflectance results are reported as percentage (%) of the total number of incident photons. 
Figure 1 illustrates the delivery/collection geometry employed in the simulations. 

Table 1. Optical parameter values used in the MC simulations 
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Fig. 1. Illumination/collection geometry employed in the simulations. 



2.2. Data analysis 

Data analysis and fitting was performed using the TableCurve2D software package (Systat) 
which utilized the Levenberg-Marquardt minimization method. TableCurve2D enabled the 
identification of the simplest and more appropriate functions to describe the data by fitting the 
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data to a library of thousands of trial functions and categorizing the results according to the 
goodness of fit. The exact functional form of Eqs. (3)-(8) was determined in this way. 

3. Results 

3.1. General model 

Diffuse reflectance calculated using the MCML code for various detector diameters is shown 
in Fig. 2 as a function of the reduced scattering coefficient. Note that for smaller detector 
diameters, reflectance exhibits an approximately linear dependence on the reduced scattering 
coefficient and this holds for all values of the absorption coefficient. The solid lines in Fig. 2 
represent a fit to the MC data using the model described by Eq. (1), 

# = — L (1) 

In Eq. (1), R is the reflectance while fi a and fi' are the absorption and reduced scattering 

coefficients, respectively. This is a simple expression that can be derived from basic physics 
[14] for the description of total reflectance (i.e. infinite detector radius). Here, the parameters 
k x and k 2 depend on the absorption coefficient and detector radius such that Eq. (1) can be 
used to describe the reflectance collected by an optical detector with arbitrary radius. These 
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Fig. 2. Reflectance as a function of the reduced scattering coefficient (MC data) for various 
detector radii, r. and various values of the absorption coefficient, (a) /u a = 0.001 mirT 1 , (b) fi a 

= 0.01 mnT 1 , (c) /J a =0.1 mnT 1 , and (d) /J a = 1 mnT 1 . Solid lines represent a fit to the model 
described by Eq. (1). 
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Probe Radius (irm) Probe Radius (mm) 

Fig. 3. k\ and k 2 parameter values obtained from the fitting results shown in Fig. 2, using Eq. 
(1). The solid lines represent spline interpolation between the data points. 
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Fig. 4. Reflectance as function of the absorption coefficient for detectors with various radii. 
The solid lines represent fits to the MC data using Eq. (2) ( fx\ = 1.8 mnT 1 for the data shown 
here). 



parameters can be determined by fitting the MC data to Eq. (1) (Fig. 2). Figure 3 shows the 
resulting values of k x and k 2 as a function of detector radius for different values of the 
absorption coefficient. Overall, the dependence of both parameters on the detector radius and 
the absorption coefficient is quite smooth and once partially determined, both k t and k 2 can be 
fully determined using spline interpolation and a lookup table technique [15] to evaluate the 
reflectance collected using an optical detector of arbitrary radius over the entire range of 
optical parameter values relevant for biological tissue. 

A more straightforward modeling approach can be identified by focusing on the 
dependence of reflectance on absorption. Figure 4 shows the reflectance for detectors with 
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various diameters as a function of the absorption coefficient with the solid lines representing 
MC data fits to Eq. (2), 

R = — ^— (2) 

As can be seen, Eq. (2) describes the MC data very well. All data shown in Fig. 4 are for 
ju' s =1.8 mm -1 with different values of the reduced scattering coefficient producing similar 
results. This confirms that the reflectance when using an optical detector with finite diameter, 
exhibits a simple inverse dependence on the absorption coefficient. We have also previously 
discovered this dependence experimentally using a small-diameter optical probe [14]. This 
observation simplifies modeling of reflectance considerably. What is needed in order to 
completely determine a general reflectance model based on Eq. (2) is to establish the exact 
dependence of the two parameters k l and k 2 on the reduced scattering coefficient and the 
detector radius r, i.e., k l {p[ , r) and k 2 (ju' s , r) . 

Figure 5 shows the values of these last two parameters as a function of the detector radius. 
The solid lines in Fig. 5 represent fits to Eqs. (3) and (4) for k t and k 2 , respectively: 
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Fig. 5. Dependence of the k { and k 2 parameters of Eq. (2) on the detector radius. The solid lines 
represent fits to Eqs. (3) and (4) for k\ and k 2 , respectively. 

Equations (3) and (4) include four parameters to be determined: a x ,b x , a, and b 2 . The last 
step in fully developing the model is to determine the dependence of these four parameters on 
the reduced scattering coefficient. It turns out that the dependence is again simple enough 
such that concise expressions can be found describing these parameters (Fig. 6 and Eqs. (5)- 
(8)): 

a, = (c, +d 1 1 (5) 



b^^+fxr 3 ' 2 )- 1 
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b 2 =e 2 + fX 



(6) 
(7) 
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Thus, an analytical expression for the reflectance using a detector with arbitrary radius can 
be fully determined in terms of Eq. (2) and eight empirically determined parameters which are 
summarized in Table 2. This expression combines the simplicity of a concise analytical 




Fig. 6. Dependence of the parameters a x ,b x , a 2 and b 2 of Eqs. (3) and (4) on the reduced 
scattering coefficient. The solid lines represent fits to Eqs. (5)-(8). 

expression with the accuracy of the MC method. It is important to point out that due to the 
semi-empirical nature of the model it is only valid within the optical parameter ranges 
specified in Table 1: fi a = 0.001-10.0 mm" 1 and ju' s = 0.6-3.0 mm" 1 . 



Table 2. Model parameter values 
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3.2. Small diameter detector model 



As already mentioned (Fig. 2 above) the dependence of reflectance on the reduced scattering 
coefficient is approximately linear for small detector diameters. This effect has been 
previously observed experimentally both by us and other researchers [3,4,14,16] and 
simplifies greatly modeling of diffuse reflectance when using small diameter detectors. To 
model reflectance from such a small diameter detector, an even simpler analytical expression 
can be found described by Eq. (9): 
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Figure 7 shows that Eq. (9) can describe the reflectance when using a detector with 0.5 
mm radius remarkably well. In contrast to the general model developed in the previous 
section, the parameters k\ and k 2 have fixed constant values and they are independent of the 
reduced scattering coefficient (k t = 6.1 and k 2 = 1.55 in this case). In addition, Fig. 7(a) 
demonstrates again here the inverse dependence of reflectance on the absorption 




Fig. 7. Diffuse reflectance scored using an optical detector with a radius of 0.5 mm showing the 
dependence of reflectance on (a) the absorption coefficient (reflectance shown is normalized to 
unity at the smallest value of the absoiption coefficient), and (b) the reduced scattering 
coefficient. Solid lines represent fits to Eq. (9). 

coefficient, which was also pointed out in the previous section (Fig. 4). It is interesting to note 
that we have experimentally observed previously that Eq. (9) describes very well the 
reflectance measured with a standard six-around-one 200 micron optical fiber probe [14]. 

Another interesting observation can be pointed out here regarding the inverse dependence 
of reflectance on the absorption coefficient. It can be seen from Figs. 4 and 7(a) that this 
dependence (expressed by Eqs. (1), (2), and (9)) describes the reflectance remarkably well for 
detector radii smaller than 5 mm but starts to deviate from the MC data for larger detector 
diameters. In this large detector diameter case, and also for a detector with infinite diameter 
(total reflectance collected) the MC data can be better described by Eq. (10) 

/ j = jt ^W<) (10) 

Figure 8 shows an example of the excellent description that Eq. (10) provides for the total 
reflectance (r = co) compared to Eq. (2), which does not perform that well for large detector 
diameters and total reflectance. Equation (10) is actually very similar to that obtained 
analytically for the total diffuse reflectance using the path integral method (see review in [14] 
and references therein). The parameter values used to obtain the fits shown in Fig. 8 were k x = 
82.66 and k 2 = 15.65 in the case of Eq. (2) and k x = 94.96 and k 2 = 3.31 in the case of Eq. (10). 
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Fig. 8. Reflectance for a detector with infinite radius and best fits using Eqs. (2) and (10). 
Equation (10) gives a better fit compared to Eq. (2), and this is also true for large detector 

diameters ( « =1.8 mirf' for the data shown here). 
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Fig. 9. Effects of detector numerical aperture as a function of the absorption coefficient for 
three different values of the reduced scattering coefficient: (a) ft' = 0.6 mm" 1 , (b) ft' =1.8 

mm -1 , and (c) ft' = 3.0 mm -1 . The efficiency parameter shown in the vertical axis represents 

the percent ratio of reflectance scored in the 0-15 deg exit angle range over the total reflectance 
scored over the entire 0-90 deg exit angle range. Solid lines represent spline interpolation 
between data points. 



#152212 - $15.00 USD Received 2 Aug 201 1; revised 6 Nov 201 1; accepted 8 Nov 201 1; published 9 Nov 201 1 
(C) 201 1 OSA 1 December 201 1 / Vol. 2, No. 12 / BIOMEDICAL OPTICS EXPRESS 3292 



3.3. Angular dependence effects 

All analysis and modeling work presented in the previous sections is based on the reflectance 
integrated over all exit angles. Some optical detection schemes such as those employing 
optical fibers typically collect diffuse reflectance within a limited exit angle range depending 
on the actual numerical aperture of the optical fibers used in probe construction as well as the 
tissue refractive index. Thus, maximum collection angles may lie typically in the 10-20 
degrees range, measured from the perpendicular to the tissue surface, while other detectors 
such as integrating spheres are theoretically capable of collecting all light exiting up to angles 
of 90 degrees. Most of the published diffuse reflectance models neglect exit angle effects [3- 
5] and in this section, we set out to evaluate such effects using our MC results. 

Figure 9 shows the differences between a fiber optical probe detector and an integrating 
sphere detector. The vertical axis "efficiency" shown in Fig. 9 is the ratio of light collected 
within the 0-15 degree exit angle range versus the full 0-90 degree range. As can be seen, this 
efficiency is approximately constant for small values of the absorption coefficient but starts 
dropping when absorption becomes comparable to the reduced scattering coefficient. But even 
in the high absorption case the spectral distortions introduced are not higher than about 15% 
and they are more prominent for smaller values of the reduced scattering coefficient, as 
expected. 

4. Discussion 

The problem of diffuse reflectance from a semi-infinite turbid medium has been addressed by 
many researchers in the field of biomedical optics. However, to the knowledge of the authors, 
there is no detailed treatment available based on MC calculations for a detector scheme with 
arbitrary diameter. The modeling approach presented here is characterized by its accuracy (an 
attribute of the MC method), its generality, and its simplicity. In addition, it goes beyond the 
restrictions of the MC method by providing a simple semi-empirical analytical expression for 
the reflectance. 

The modeling geometry employed in this work is potentially useful for modeling realistic 
complex optical probe geometries. A similar approach has been successfully employed in the 
past [17] to model reflectance collected by fiber optic probe (using an "effective" maximum 
radius that best models the probe). However, that work was based on diffusion theory [5]. 
Compared to that work, the present work presents significant advantages mainly due to the 
fact that it is not limited by the inaccuracies and applicability limitations of the diffusion 
approximation. 

It is interesting to compare the present work with another MC -based modeling technique 
[18]. This scaling MC method offers significantly faster computational times without 
sacrificing the modeling accuracy associated with standard MC methods. The advantage of 
our approach compared to such a method is the use of an analytical expression, albeit semi- 
empirical in nature, while retaining the accuracy associated with the MC method. The 
disadvantages of our approach lie in the restrictions associated with the specific 
illumination/collection geometry and the semi-infinite tissue approximation. However, these 
disadvantages do not seem to impose serious limitations because the modeling approach 
presented here can be potentially extended to more general delivery/collection geometries as 
well as multilayer tissue geometries by following a similar methodology in model 
development. We hope that we will be able to report on such generalizations of the present 
modeling approach in the future. 

The specific delivery/collection geometry of the present model (i.e. an infinitely thin 
delivery beam and idealized detector with fixed radius) may at first appear far from realistic 
delivery/collection schemes. Nevertheless, we must point out that it may be quite useful in 
many practical situations. As already mentioned above, a similar approach has been 
successfully applied using a fiber-optic probe [17] and other similar approaches have also 
been reported [14,19,20]. In addition, the present approach may be useful in modeling 
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integrating sphere detectors and/or delivery/collection geometries with free space illumination 
and a camera for reflectance collection. 

Regarding the similar previous published models, it is worth mentioning those that employ 
a modeling approach similar to that described by Eq. (9) ([14,21-25] by us and [19,20] by 
other researchers) to describe reflectance collected using optical fiber probes. In particular, the 
simple inverse dependence of reflectance on the absorption coefficient (illustrated in Figs. 4 
and 7(a)) is worth mentioning because it greatly simplifies modeling and, to the authors' 
knowledge, it has not been previously reported by other researchers. This, in combination 
with the simple linear dependence of reflectance on the reduced scattering coefficient (which 
has been reported and confirmed by other researchers previously) makes for the very simple 
expression described by Eq. (9). Equations (1), (2), and (9) may indeed seem very simple and, 
in fact, they may even look too simplistic to some researchers. Nevertheless, they appear to be 
quite promising and effective and we would like to urge researchers in the field of biomedical 
optics to investigate their utility further. At least two other publications [19,20] have already 
reported on the successful application of a reflectance model similar to that described by Eq. 
(9). 

Equation (1) relies on a look-up table technique which has already been suggested by 
other researchers as a convenient modeling tool in tissue reflectance [15]. In addition, it is 
possible to substitute the lookup table technique with an empirical fit of the data as is 
illustrated for Eq. (2). The problem in such a case is to identify the appropriate analytical 
expression that best describes the data. We have shown here that it is possible to develop such 
a semi-empirical model based on Eq. (2) using this approach, as we have also done previously 
for a two-layer reflectance model [25]. In both cases, the complexity of the model is 
significantly lower and the ease of model application and accuracy are potentially greater 
compared to other available reflectance models such as those based on diffusion theory [5,17]. 

Investigation of the angular dependence of diffuse reflectance revealed that the finite 
numerical aperture of a fiber optic probe does not introduce any appreciable spectral 
distortions as long as absorption is lower that scattering. Even in the limit of high absorption 
(which is sometimes present in biological tissues) the error introduced is no larger than about 
15%. These results are consistent with the practice of many researchers who neglect exit angle 
effects in their modeling by using models that are directly based on the total reflectance at all 
exit angles [3-5]. When accurate measurements are required though, and absorption is high, 
the angular dependence of reflectance must be taken into account. 

Another important point is that the same modeling approach described here using Eqs. (1), 
(2), and (9) can be applied to experimental data on tissue phantoms. In this way, a more 
accurate, reliable, and realistic calibration of the model can be achieved by either tuning the 
actual model parameter values or even coming up with new model analytical expressions such 
as those described by Eqs. (1), (2), (9), and (10). We have already employed this modeling 
approach using tissue phantoms in the development of both one and two-layer tissue 
reflectance models [14,21,25] and we hope that this approach will serve as the basis for more 
accurate and robust modeling of tissue reflectance in the future. 

5. Conclusions 

This work presented a simple semi-empirical tissue reflectance model useful for an optical 
detector with arbitrary diameter. The development of the model is based on MC simulations 
results and thus it is not restricted by the limitations of the diffusion approximation. In 
addition, this modeling approach utilizes two simple characteristics of the reflectance: (1) an 
inverse dependence on the absorption coefficient, and (2) a linear dependence on the reduced 
scattering coefficient (for small diameter detectors). We expect that the modeling approach 
presented here will facilitate further tissue analysis and investigations using reflectance 
spectroscopy and that it will be useful to researchers in the field of biomedical optics. 
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